clear; close all; addpath('../utilities');

para = [0.011811964656	0.244090504804	0.034225636514	0.035428703735	2.565573269391	0.165210276777	0.298220369064];

main_para

options = optimoptions('fsolve','Display','iter','TolFun',1e-12);
beta    = fsolve(@solve_beta,0.99,options,gamma,psi,mu_C,sig_C,P);

main_SDF

lamS    = ones(nS,nK);
lamD    = ones(nS,nK);
k       = linspace(k_min,k_max,nK)';
k_bar   = ones(nS,1)*5;
k_iss   = ones(nS,1)*10;

[kgrid,mgrid]   = meshgrid(k,(1:nS)');
kgrid           = kgrid(:);
mgrid           = mgrid(:);

Q_big   = repmat(Q,nK,1);
Rf1_mat = repmat(Rf(:,1),1,nK);

main_VFI

main_CDS

main_CDS_P

money       = ones(nS,nK); % moneyness levels

main_options

IV          = NaN(nS,nK,2);
IV(:,:,1)   = BSIV;

money       = exp(-BSIV/sqrt(12)); % moneyness levels

main_options

IV(:,:,2)   = BSIV;
IV_ATM      = IV(:,:,1)';
IV_ATM      = permute(IV_ATM,[2,1]);
IV_Skew     = IV(:,:,2)'-IV(:,:,1)';
IV_Skew     = permute(IV_Skew,[2,1]);

main_Ret

main_simu

save model_simulation table_5 G_m
